tile=readLAS("../tests/data/NEON_D03_OSBS_DP1_405000_3276000_classified_point_cloud.laz")
tile<-tile %>% lasfilter(Z < 100)
tile<-tile %>% lasfilter(Z > 0)
plot(tile, color="Z", colorPalette = heat.colors(50), trim=0.95)
#remove
summary(tile)
#table(tile@data$Classification)
#ggplot(tile@data,aes(x=as.factor(Classification),y=Z)) + geom_boxplot()
#add a point index
tile@data$PointID<-1:nrow(tile@data)
#Lastrees updates the cloud inplace, save data in seperate object
ptlist<-list()
#ground model
ground_model(las)
plot(li2012,color="treeID",colorPalette = sample(pastel.colors(500)), size = 1)
#canopy model
chm=canopy_model(las)
#Toy file for testing
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
tile = readLAS(LASfile, select = "xyz", filter = "-drop_z_below 0")
col = pastel.colors(200)
#add a point index
tile@data$PointID<-1:nrow(tile@data)
#Lastrees updates the cloud inplace, save data in seperate object
ptlist<-list()
ground_model(tile)
chm=canopy_model(tile)
plot(tile,color="Classification")
You must enable Javascript to view this page properly.
system.time(watershed<-segment_ITC(tile,algorithm = "watershed",chm=chm))
## user system elapsed
## 5.438 0.098 19.046
ptlist[["watershed"]]<-watershed %>% lasfilter(!is.na(treeID)) %>% .@data
plot(watershed,color="treeID",colorPalette = sample(pastel.colors(500)), size = 1)
system.time(li2012<-segment_ITC(tile,algorithm = "li2012",chm=chm))
## user system elapsed
## 0.284 0.025 1.217
ptlist[["li2012"]]<-li2012 %>% lasfilter(!is.na(treeID)) %>% .@data
plot(li2012,color="treeID",colorPalette = sample(pastel.colors(500)), size = 1)
system.time(dalponte2016<-segment_ITC(tile,algorithm = "dalponte2016",chm=chm))
## user system elapsed
## 4.301 0.051 17.771
ptlist[["dalponte2016"]]<-dalponte2016 %>% lasfilter(!is.na(treeID)) %>% .@data
system.time(silva2016<-segment_ITC(tile,algorithm = "silva2016",chm=chm))
## user system elapsed
## 4.268 0.054 17.102
ptlist[["silva2016"]]<-silva2016 %>% lasfilter(!is.na(treeID)) %>% .@data
library(pander)
df<-data.frame(Algorthm=c("watershed","li2012","Dalponte2016","Silva2016"),Total_Trees=c( max(ptlist[["watershed"]]$treeID,na.rm=T),max(ptlist[["li2012"]]$treeID,na.rm=T), max(ptlist[["dalponte2016"]]$treeID,na.rm=T), max(ptlist[["silva2016"]]$treeID,na.rm=T)))
pandoc.table(df,style="rmarkdown")
##
##
## | Algorthm | Total_Trees |
## |:------------:|:-----------:|
## | watershed | 76 |
## | li2012 | 378 |
## | Dalponte2016 | 111 |
## | Silva2016 | 111 |
pdf<-melt(ptlist,id.vars=colnames(ptlist[["silva2016"]]))
head(pdf)
## X Y Z PointID Classification Zref treeID L1
## 1 481352.0 3813013 0.000 23 2 0.33 2 watershed
## 2 481352.0 3813013 0.000 24 2 0.60 2 watershed
## 3 481352.3 3813015 14.289 25 0 15.18 2 watershed
## 4 481352.3 3813015 14.835 26 0 15.49 2 watershed
## 5 481352.2 3813014 14.595 27 0 15.02 2 watershed
## 6 481352.2 3813015 14.409 28 0 14.97 2 watershed
Since the input of segmentation methods differ among algorithms, there may be different thresholds for what constitutes a tree. Let’s remove points that only appear in less than 3 methods.
#How many methods does each point appear in.
table(table(pdf$PointID))
##
## 1 2 3 4
## 5342 2070 2976 35239
points_to_remove<-pdf %>% group_by(PointID) %>% summarize(n=n()) %>% arrange(n) %>% filter(n < 3) %>% .$PointID
pdf<-pdf %>% filter(!PointID %in% points_to_remove)
Did that improve the similarity of # of tree objects?
pdf %>% group_by(L1) %>% summarize(Total_Trees=length(unique(treeID)))
## # A tibble: 4 x 2
## L1 Total_Trees
## <chr> <int>
## 1 dalponte2016 111
## 2 li2012 356
## 3 silva2016 111
## 4 watershed 76
#Create point ID lookup table
res<-dcast(pdf,PointID~L1,value.var = "treeID")
idframe<-res %>% add_rownames() %>% select(rowname,PointID)
head(res<-res %>% select(-PointID))
## dalponte2016 li2012 silva2016 watershed
## 1 110 35 5 2
## 2 110 35 17 2
## 3 111 35 5 2
## 4 111 35 5 2
## 5 110 35 5 2
## 6 111 35 5 2
system.time(consensus<-majority_voting(res, is.relabelled = FALSE))
## user system elapsed
## 6.648 0.105 20.573
#reassign to pointID
consensus_frame<-data.frame(rowname=rownames(res),consensus=consensus) %>% inner_join(idframe) %>% select(PointID,consensus)
#merge with the original tile
tile@data<-merge(tile@data,consensus_frame,all=T)
plot(tile,color="consensus",colorPalette = sample(pastel.colors(500)), size = 1)
You must enable Javascript to view this page properly.
but let’s just look at some individual trees. The question is, can we classify them based on the lidar cloud? Do we need to drape on the RGB data.
ind_trees= split(tile@data, tile@data$treeID)
ind_trees = lapply(ind_trees, LAS, header = tile@header)
plot(ind_trees[[20]],color="treeID",bg="grey90")
You must enable Javascript to view this page properly.
plot(ind_trees[[30]],color="treeID",bg="grey90")
You must enable Javascript to view this page properly.
plot(ind_trees[[40]],color="treeID",bg="grey90")
You must enable Javascript to view this page properly.
Clearly, more work to be done. Especially in reference to the ground model.
#save.image("ITCs.RData")